#include "libmri.h"
#include "auxfns.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <omp.h>

#include <iostream>
using std::cout;
using std::endl;
using std::cerr;

extern void pocs2d_l1spirit_omp (cplxf *data, int *dsize, cplxf *k2d,
                                 l1spirit_config const& config);

extern void pocs2d_l1spirit_cuda (cplxf *data, int *dsize, cplxf *k2d,
                                  l1spirit_config const& config);

void l1spirit (cplxf *data, int *dsize, l1spirit_config const& config)
{
	cplxf *calib;
	int csize[4];

	cplxf *k3d, *k2d;
	int *ksize = config.ksize;

/***************
 * Scale the data to allow re-use of tuned algorithm parameters
 */
	float scale_factor;
	rescale_data (data, dsize, scale_factor, config);

	l1spirit_debug(config,data,0,data,4,dsize[0],dsize[1],dsize[2],dsize[3]);
	l1spirit_debug(config,scale,0,&scale_factor,1,1);

/***************
 * Parallel Imaging Calibration:
 * - Crop the calibration region into the calib array
 * - Solve the calibration equations via the appropriate calib3d_subroutine
 */
	crop_calib_data (calib, csize, data, dsize, config.calib_n_readouts);

	l1spirit_debug (config,calib,0,calib,4,csize[0],csize[1],csize[2],csize[3]);

	double ts = timestamp();
	switch (config.calib_algorithm)
	{
		case l1spirit_config::lanczos: 
		case l1spirit_config::culanczos:
		{
			int nch = dsize[3];
			int n = ksize[0]*ksize[1]*ksize[2]*nch;
			int k = n;
			int nv = -1;
			float tau = config.calib_threshold;
			cplxf *V;

			if (config.calib_algorithm == l1spirit_config::culanczos){
				cout << "calib3d_culanczos()" << endl;
				calib3d_culanczos (V, nv, calib, csize, ksize, k, tau);
			}
			else{
				cout << "calib3d_lanczos()" << endl;
				calib3d_lanczos (V, nv, calib, csize, ksize, k, tau);
			}

			int kksz[3] = {2*ksize[0]-1, 2*ksize[1]-1, 2*ksize[2]-1};
			k3d = new cplxf [kksz[0]*kksz[1]*kksz[2]*nch*nch];

			double tskr = timestamp();

			if (config.calib_algorithm == l1spirit_config::culanczos){
				cout << "cukreduce3d()" << endl;
				cukreduce3d (k3d, V, ksize, nch, nv);
			}
			else{
				cout << "kreduce3d()" << endl;
				kreduce3d (k3d, V, ksize, nch, nv);
			}
			cout  << "kreduce: " << (timestamp()-tskr) << " s" << endl;

			for (int i = 0; i < 3; ++i)
				ksize[i] = kksz[i];

			delete [] V;
			break;
		}

		case l1spirit_config::cheev: {
			cout << "calib3d_cheev()" << endl;		

			complex *V = 0;
			int nv = -1;
			int nch = dsize[3];

			float tau = config.calib_threshold;

			calib3d_cheev (V, nv, (complex*) calib, csize, ksize, tau);

			int kksz[3] = {2*ksize[0]-1, 2*ksize[1]-1, 2*ksize[2]-1};
			k3d = new cplxf [kksz[0]*kksz[1]*kksz[2]*nch*nch];

			kreduce3d (k3d, (cplxf*) V, ksize, nch, nv);

			for (int i = 0; i < 3; ++i)
				ksize[i] = kksz[i];

			delete [] V;

			break;
		}

		case l1spirit_config::posv:
			cout << "calib3d_posv()" << endl;
			calib3d_posv (k3d, ksize, calib, csize, config.calib_tychonov);
			break;


	  default: default_calib_algorithm:

		case l1spirit_config::chol:{
			cout << "calib3d_chol()" << endl;
			calib3d_chol (k3d, ksize, calib, csize, config.calib_tychonov);
			break;
		}

	};
	l1spirit_debug(config,calib_k3d,0,k3d,5,ksize[0],ksize[1],ksize[2],csize[3],csize[3]);

	ifft_kern_readout (k2d, k3d, ksize, dsize);

	l1spirit_debug(config,calib_k2d,0,k2d,5,dsize[0],ksize[1],ksize[2],csize[3],csize[3]);
	cout << "Calibration: " << (timestamp()-ts) << endl;



/**************
 * Compressive Sensing l1 Minimization
 * 1 Inverse Fourier Transform along the readout dimension to decouple the
 *   2D problems. Do this both for the data itself and the calibration kernels
 * 2 Solve the 2D problems via POCS (in pocs2d_l1spirit)
 * 3 Fourier Transform the readout dimension to return the data to  k-space
 */

	ifft_readout (data, dsize);

	l1spirit_debug(config,data_roifft,0,data,4,dsize[0],dsize[1],dsize[2],dsize[3]);

	if (config.l1min_on_gpu)
		pocs2d_l1spirit_cuda (data, dsize, k2d, config);
	else
		pocs2d_l1spirit_omp (data, dsize, k2d, config);

	l1spirit_debug(config,data_soln,0,data,4,dsize[0],dsize[1],dsize[2],dsize[3]);

	fft_readout (data, dsize);

	l1spirit_debug(config,data_rofft,0,data,4,dsize[0],dsize[1],dsize[2],dsize[3]);

/***************
 * Return the output to the same scale as the input
 */
	unscale_data (data, dsize, scale_factor);

	delete [] calib;
	delete [] k3d;
	delete [] k2d;
}

